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Ozone and particulate matter PM2.5 are co-pollutants that have 
long been associated with increased public health risks. Information 
on concentration levels for both pollutants come from two sources: 
monitoring sites and output from complex numerical models that pro- 
duce concentration surfaces over large spatial regions. In this paper, 
we offer a fully-model based approach for fusing these two sources 
of information for the pair of co-pollutants which is computationally 
feasible over large spatial regions and long periods of time. Due to 
the association between concentration levels of the two environmen- 
tal contaminants, it is expected that information regarding one will 
help to improve prediction of the other. Misalignment is an obvious 
issue since the monitoring networks for the two contaminants only 
partly intersect and because the collection rate for PM2.5 is typically 
less frequent than that for ozone. 

Extending previous work in Berrocal et al. (2010), we introduce a 
bivariate downscaler that provides a flexible class of bivariate space- 
time assimilation models. We discuss computational issues for model 
fitting and analyze a dataset for ozone and PM2.5 for the ozone sea- 
son during year 2002. We show a modest improvement in predictive 
performance, not surprising in a setting where we can anticipate only 
a small gain. 

1. Introduction. Ozone and particulate matter, PM2. 5, have long been 
associated with increased public health risks, e.g., of respiratory diseases 
[4, 12, 45], cardiovascular diseases [4, 12], and mortality and morbidity in 
general [13, 47]. To set air quality standards the EPA utilizes information 
from monitoring networks and from air quality computer models. 

The two sources of information are valuable in different ways. The ob- 
servations reported by monitoring networks, though sparsely collected and, 
sometimes, affected by missingness, provide direct measurements of the true 
value up to measurement error. The output from air quality models esti- 
mates pollutant concentrations over large spatial domains at the grid cell 
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level. The estimates are viewed as averages over these cells and do not con- 
tain any missing data, but are uncalibrated. It is desirable to combine both 
sources of information but to do so, it is necessary to solve the inherent 
"change of support" problem of the data [2, 9, 22], that is, the spatial mis- 
alignment between the observational data and the numerical model output. 
Addressing the difference in spatial resolution between the two sources of 
data will also allow to evaluate and calibrate the numerical model output. 

We are not interested in calibration of the numerical model in the sense 
of [29]. Rather our focus is on resolving the difference in spatial scale be- 
tween the numerical model output and the observational data, while bias- 
correcting the predictions generated by a numerical model. In this regard, 
several methods have been proposed to assess numerical models and address 
the issue of "incommensurability" between grid model averages and point 
measurements [50]. Such effort is also common in the context of data as- 
similation within the atmospheric science literature, where the goal is to 
combine observational data on the current state of the atmosphere with a 
short-range forecast produced by a numerical weather prediction model to 
yield initial conditions for a numerical atmospheric model. The approaches 
tend to be algorithmic and ad hoc, only occasionally model-based, and do 
not necessarily address downscaling. Fuller discussion can be found in [10] 
and [28]. In the context of air quality, [36] propose to evaluate the hourly 
ozone predictions generated by a geophysical model on the grid-cell scale 
by using the observations taken at monitoring sites, and predicting hourly 
level of ozone at the model grid cells by averaging the predictions at M 
regularly-spaced sites taken within each grid cell (see, as well, [15]). A dif- 
ferent strategy has been proposed by [27] who ignore the difference in spatial 
resolution between model output and observations, rather suggesting to eval- 
uate a numerical model by looking for differences between the model output 
and the observations in terms of variograms and correlograms. 

[16] develop a "fusion" model, expressing the numerical model output as 
the integral over a grid cell (scaled by the area of the cell) of a latent point- 
level process. Their main goal is to recover the true unobserved process, 
but, as a by-product, they obtain estimates of the bias of the numerical 
model output, thus allowing evaluation and calibration of the numerical 
model. The work is an instance of Bayesian melding [39], and has gained 
considerable popularity [14, 46]. See also recent work in this spirit applied 
to sulfate aerosol [49] and ammonium [11]. The Bayesian melding approach 
of [16], though popular, suffers from several limitations, as pointed out by 
[33] and [3]. Additionally, it is only computationally feasible as a spatial 
model. A spatio-temporal extension, proposed by [35], specifies the latent 
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process at the grid cell level, leading to comparison between the numerical 
model output and the observational data at the grid cell scale. 

Operating at the point level, [23] propose to correct the model output by 
downscaling the predictions from grid cells to points and comparing them 
with the observations. For each site, a two-step linear regression model that 
relates the observed hourly ozone level at a site with the numerical model 
output at the grid cell in which the site lies is proposed. Building upon 
the work of [23] , [32] suggest interpolating the intercept and the coefficient 
of the numerical model output, estimated at each site, via kriging, thus 
developing a "spatio-temporal" model that will allow corrected prediction 
of ozone level also at unmonitored sites. However, this ad-hoc approach does 
not provide estimates of the uncertainties associated with the predictions. 
A formal spatio-temporal model, that builds upon the downscaling idea of 
[23], has been proposed as a univariate downscaler by [3] with a version 
implemented and described in [42] and in [18]. 

The contribution of this paper is to extend the spatio-temporal downscaler 
from the univariate setting to a bivariate setting in order to exploit the cor- 
relation both in the observed concentration levels of ozone and PM2.5 and 
in the concentration levels of ozone and PM2.5, provided by the numerical 
model. Working in a bivariate setting will prove particularly advantageous 
for predicting particulate matter. The sampling frequency of the PM2.5 mon- 
itors - most of the monitors measure PM2.5 concentration every 3 days - 
yields challenging interpolation of the monitoring data to the entire spatial 
domain as well as difficult evaluation of the predictions generated by the nu- 
merical model. Rather, in environmental health studies, researchers always 
use monitoring data to characterize exposure to fine particulate matter, 
dealing with the PM2.5 sampling frequency issue by aggregating the mon- 
itoring data over time and/or over space. Similarly, most of the research 
effort on PM2.5 has been devoted to developing spatial and spatio-temporal 
models to predict PM2.5 based on monitoring data, meteorological covari- 
ates or observed concentrations for another pollutant [5, 30, 31, 40, 43, 48]. 
A fusion model to combine monitoring data for PM2.5 with satellite aerosol 
optical depth (AOD) data has been proposed by [37]: however, the analysis 
was conducted on monthly average concentration rather than daily and it 
revealed that AOD provides no improvement in predicting fine particulate 
matter. 

In this paper, we present a general bivariate spatio-temporal downscaler 
developed for the numerical model outputs of an air quality model; we illus- 
trate the methodology with regard to prediction of concentration of ozone 
and fine particulate matter. The model, which to our knowledge, is the 
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first bivariate fusion/downscaler model, extends the univariate downscaler 
of [3], and, following [3], regresses the bivariate vector of observed ozone 
and PM2.5 concentration on the numerical model outputs for both pollu- 
tants using spatially- varying coefficients [19, 44], modeled as a correlated 
six-dimensional Gaussian process. The model specification is general and 
flexible, while offering feasible computation for a fully model-based fusion 
across space and time. We explore several different special cases, correspond- 
ing to different assumptions on the correlation structure of the two pollu- 
tants. An important feature of our model is that it allows us to handle not 
only the spatial misalignment between monitoring data and model output, 
but also it allows us to accommodate the spatial and temporal misalign- 
ment between the ozone and PM2.5 monitoring data. Finally, exploiting the 
correlation, not only between the model output and the corresponding mon- 
itoring data for a pollutant, but also between pollutants, we can jointly 
predict ozone and PM2.5 at any site in the spatial domain and provide a 
measure of the uncertainties associated with such predictions. Additionally, 
at sites where one pollutant is measured but the other is not, we can predict 
the level of the latter. 

The paper is organized as follows. In Section 2, we present the moni- 
toring data and the numerical model output, in Section 3 we describe the 
bivariate downscaler model, by first introducing the downscaler in the uni- 
variate setting, and then extending the model to the bivariate case, first in a 
purely spatial setting and then in a spatio-temporal setting. Details on how 
to handle the spatio-temporal misalignment in the monitoring data are also 
discussed in Section 3. Section 4 presents details on the model fitting, while 
Section 5 displays results of our analysis. Finally, we conclude in Section 
6 with a discussion and evaluation of our method and with indications for 
future work. 



2. Data. Particulate matter (PM2.5) and ozone are two of the "criteria 
pollutants" that the Environmental Protection Agency (EPA) is required 
to monitor by the Clean Air Act. The first is a mixture of solid and liquid 
particles, emitted in the atmosphere either directly from a source or as result 
of complicated reactions of chemicals, while the second is a gas made of three 
atoms of oxygen that is created by a chemical reaction between oxides of 
nitrogen and volatile organic compounds in the presence of sunlight. 

The EPA tracks both criteria pollutants using both measurements taken 
at monitoring sites and estimates of air pollutants concentration produced by 
the Models-3/Community Mesoscale Air Quality (CMAQ) model [6] (epa.gov/asmdnerl/CMAQ). 
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The latter is a deterministic numerical model that predicts concentrations for 
various pollutants by integrating three major components: an atmospheric 
component accounting for the atmosphere and its states and motions, an 
emission component accounting for the emissions injected in the atmosphere, 
and a chemical component accounting for the reactions between the different 
gas present in the atmosphere. By simulating various chemical and physi- 
cal processes, such as horizontal and vertical advection, emission injection, 
deposition, plume chemistry effects, etc., CMAQ produces estimates of pol- 
lutants concentration at pre-determined spatial scales, e.g., 36 km, 12 km, 
and, most recently, at 4 km grid cells. 

We consider daily estimates and measurements of ozone and PM2.5 con- 
centration for the region in the South Eastern part of the United States 
shown in Figure 1 during the period June 1 - September 30, 2002, the sum- 
mer season when solar radiation and temperature are high. In turn this en- 
courages not only high concentrations of ozone but also of particulate SO4, 
the dominant component of PM2.5 in the eastern U.S as well as particulate 
ammonium nitrate. Relatively strong association between ozone and PM2.5 
concentrations is expected and, in fact, observed. Following the air quality 
standards (NAAQS; http://www.epa.gov/air/criteria.html) set by EPA for 
ozone and PM2.5, daily ozone concentration is measured as the daily maxi- 
mum 8-hour average concentration, while daily PM2.5 concentration is given 
by the 24-hour average concentration of particulate matter. 

The monitoring data used in our case study comes from 226 sites sparsely 
located in the region, and it has been obtained from the EPA Air Qual- 
ity System (AQS) repository database. Of these 226 sites, not all measure 
both pollutants: 71 report only ozone measurements, 50 report only PM2.5, 
and 105 measure both, pinpointing the spatial misalignment in the monitor- 
ing data. In order to validate the model with out-of-sample predictions, we 
randomly split the sites into two sets: a training set used to fit the model, 
and a validation set used to assess the performance of the model. In the 
training set, 52 sites measured only ozone concentration, 39 reported only 
PM2.5 concentration, and 70 reported both. In total, the dataset used to fit 
the model comprised 14,630 daily measurements of ozone concentration and 
4,790 measurements of PM2.5 concentration. The smaller number of obser- 
vations for particulate matter is due to the sampling scheme for PM2.5: of 
the 109 PM2.5 monitoring sites, only 11 (10.1%) measure concentration of 
particulate matter every day, while 90 (82.5%) report measurements every 
three days, and 8 (7.3%) measure PM2.5 every six days. 

The difference in sampling frequency between ozone and PM2.5 suggests 
that, in addition to spatial misalignment, there is also temporal misalign- 
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ment in the monitoring data: even sites that measure both pollutants do not 
sample the two pollutants with the same temporal frequency. However, when 
concentration measurements for both pollutant are available, they display 
considerable correlation, equal to 0.69. Spatial and temporal misalignment 
is also present among the validation sites: 54 measured only ozone, 46 only 
PM2.5 sites, while 35 reported concentration for both pollutants. In total, 
the validation dataset contained 6,530 daily observations of ozone concen- 
tration and 2,559 daily observations of concentration of particulate matter. 
The location of the sites used to fit the and validate the model is shown in 
Figure 1. 

An exploratory analysis of the observed daily concentration data for ozone 
and PM2.5 revealed the need for a transformation. We modeled the daily 
ozone concentration data on the square root scale to stabilize the variance, 
while we log-transformed the PM2.5 concentration data to achieve normal- 
ity. Both transformations have been previously used; see, e.g., [41], [40], [7], 
[25] . A normal Q-Q plot of the square root of the observed daily concentra- 
tion of ozone at sites used to fit the model is shown in Figure 2 (a), while 
Figure 2 (b) presents a normal Q-Q plot of the logarithm of daily average 
concentration of PM2.5 at the 109 monitoring sites in the test set measur- 
ing particulate matter. Both plots seems to suggest a slight deviation from 
normality, however the deviation is minimal and we are comfortable using 
a normal model, in agreement with the literature. Neverthless, in Section 6 
we address distributional issues more at lengths and suggest an extension 
of our downscaling approach for non-Gaussian variables. The overall mean 
and standard deviation for the square root of ozone concentration and for 
the log of PM2.5 concentration at the test sites were, respectively, 7.41 and 
1.31 \/ppb (ppb:parts per billions), and 2.72 and 0.56 log(^tg/m 3 ). 

Estimates of the daily average and standard deviation for ozone and PM2.5 
are presented in Figure 3. Both panels in Figure 3 indicate daily variability 
in the concentration level for the pollutants, with Figure 3(a) revealing the 
seasonality of ozone during the study period (June 1 - September 20, 2002). 
In terms of variability, both ozone and PM2.5 present the largest standard 
deviation on June 25, 2002. For this day, plots of the observed ozone and 
PM2.5 concentration are shown, on the original scales, respectively, in Fig- 
ure 5(a) and Figure 6(a). (These plots are developed in conjunction with 
the data analysis of Section 5.) 

The numerical model output data consists of estimates of ozone and PM2.5 
concentration level generated by the CMAQ numerical model run at 12 km 
grid cell resolution. An example of the type of concentration surfaces yielded 
by CMAQ can be observed, respectively, in Figure 5(b) and Figure 6(b). 
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CMAQ produces estimates of a pollutant concentration in terms of average 
over a grid cell (in this case, over 10,504 12-km grid cells), while observa- 
tions are taken at points. Therefore, there is a spatial misalignment between 
observational and numerical model data. We solve this difference in spatial 
resolution between the monitoring and the numerical model data by asso- 
ciating to each site s in the spatial domain S, the CMAQ grid cell B in 
which s lies. This allows a comparison between the model output and the 
monitoring data. 

To reveal the need for calibration in the numerical model output, we have 
produced pairwise scatterplots of the square root of the observed ozone 
concentration versus the square root of the CMAQ predicted concentration 
level of ozone during the period June 1 - September 30, 2002, of the square 
root of the observed ozone concentration versus the log of the CMAQ pre- 
dicted PM2.5 concentration, and analogous plots for the log of the observed 
concentration of particulate matter. We see that the predictions by the nu- 
merical model are biased and need to be calibrated; however, they do contain 
useful information to improve prediction for both pollutants. Additionally, 
the positive and substantial correlation between the CMAQ model output 
for PM2.5 and the observed ozone concentration (r=0.62), and, similarly, 
between the CMAQ model output for ozone and the observed PM2.5 con- 
centration (r=0.69), indicate that the CMAQ output for PM2.5 might be 
useful to predict ozone and conversely. 

The type of calibration implied by the pairwise scatterplots mentioned 
above is constant across the entire spatial domain S in consideration. Em- 
pirical evidence suggests instead that the error and the bias of the numerical 
model output might not be constant in space, rather they vary from site to 
site. Figure 4 presents spatial maps of the estimates of the intercept and 
of the regression coefficients of CMAQ ozone and CMAQ PM2.5 at each of 
the sites used to fit our bivariate downscaler model. These estimates have 
been obtained by regressing at each site s, respectively, the normalized log 
of the observed PM2.5 concentration at s across time on the corresponding 
normalized CMAQ model output for ozone and PM2.5, both taken on the 
appropriate scale (square root for ozone, log for PM2.5). As the plots indi- 
cate, there is spatial variability in the estimates of these coefficients with 
differential variability across the estimates. In particular, we see less spatial 
variability in the estimates obtained for the observed log PM2.5 concentra- 
tion. Also, there is a difference between the two pollutants in the significance 
of the estimates of the coefficients: while for ozone, the estimates of the coef- 
ficient of CMAQ ozone are all significant, and about 60% of the estimates of 
the intercept are significant, in the regression for PM2.5, only 21% of the in- 
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tercepts are significantly different from zero, and about 88% of the estimates 
of the coefficient of CMAQ PM2.5 are significantly different from zero. 

3. Modeling. Here we present the model underlying our downscaling 
technique. We first briefly review the univariate downscaler, presented in [3], 
then we extend it to the bivariate setting. In doing so, we first present the 
general bivariate downscaler in a purely spatial setting, then we examine 
extension to accommodate data collected over time. 

3.1. Univariate downscaler. Let Y(s) be the observed data at s and let 
x{B) be the numerical model output at grid cell B. Then, x(B) is inter- 
preted as the average level of the variable under consideration over B. The 
spatial downscaler addresses the difference in spatial resolution between the 
observational data and the numerical model output, i.e. the "change of sup- 
port" problem [2, 9], by associating to each point s the grid cell B in which s 
lies. Then, it relates the observational data and the numerical model output 
as follows: 

(1) Y(s) = /3 (s) + Mb)x(B) + e(s) e(s) # N(0, r 2 ) 
where 

/3 (s) = /?o + /3o(s) 

(2) ^(s) = /3x + A(s) 

In (2), /3q and /3i denote, respectively, the overall additive and multiplicative 
bias of the numerical model, and /?o( s ) and /3i(s) are local adjustments to 
the overall bias terms. Finally, e(s) is a white noise process with nugget 
variance r 2 . Evidently, a constant-mean downscaler is a special case of (2) 
obtained when the local adjustments, /3q(s) and /3i(s), are set to zero. 

Anticipating association between intercept and slope, at the second hier- 
archical level of the model, we set the two spatially- varying coefficients /3o( s ) 
and /3i(s) to be a bivariate correlated mean-zero Gaussian processes. Using 
the method of coregionalization [19, 44, 51], we assume that there exist two 
latent mean-zero, unit- variance independent Gaussian processes u>o(s) and 
wi(s) such that Cov(wj(s),Wj(s')) = exp(— <fij\\s — s'||), where, for j = 0, 1, 
4>j is the spatial decay parameter of the Gaussian process Wj(s), ||s — s'|| is 
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the great circle distance between s and s' 1 , and 




) 



The coregionalization matrix A in (3) determines the correlation between 
the two spatially- varying coefficients /3q(s) and /3i(s) and can be assumed, 
without loss of generality, to be lower-triangular. Note that T = A A' is 
the constant local covariance matrix. The specification of the model is then 
completed with the priors for: the overall bias terms, /3o and /3i, the nugget 
variance r 2 , the three non-null elements of the coregionalization matrix A, 
and the decay parameters 4>o and <j>\. Other choices for obtaining dependent 
Gaussian processes could be entertained, including moving averages [26], 
convolution of covariance functions [34], and particular parametric choices. 
In the last case, some possibilities include the cross-covariance functions 
introduced by [1] and [20]. 

3.2. The bivariate downscaler: static version. Turning to the bivariate 
case, we describe the model in terms of square root ozone and log PM2.5 
concentrations, recognizing that the model could be applied to any suitable 
pair of point referenced variables. Let ii(s) and ^(s) denote, respectively, 
the square root and the logarithm of the observed ozone and PM2.5 con- 
centration at s, and let x\{B) and X2(B) indicate, respectively, the square 
root and the logarithm of the numerical model output for ozone and PM2.5 
at grid cell B. Then, Y(s) = (Y"i(s), Y 2 (s)) is the bivariate random vec- 
tor with the observed concentrations for the two pollutants at s, while 
x(i?) = (xi(B), x 2 (B)) is the bivariate vector with the two numerical model 
outputs. Building upon (1), our bivariate downscaler model relates the mon- 
itoring data at s and the numerical model outputs at grid cell B, with s lying 
in B, as follows: 



where ei(s) and e2(s) are two independent white noise processes with nugget 
variances, respectively, r 2 and r| • 

1 When considering large spatial domains, it is preferable to use three-dimensional Eu- 
clidean distance as an argument for the exponential covariance function rather than great- 
circle distance, which might yield a covariance matrix which is not positive definite. For 
the domain we study, the two metrics almost coincide and yield very similar results. 



(4) 



Y 2 (s) 



$io(b) + /9n(s)a;i(B) + Pu(b)x 2 (B) + ei(s) 
/3 20 (s) + ki(s)x!(B) + p 22 (s)x 2 (B) + e 2 (s) 
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As in the univariate case, we adopt a random effect notation, decomposing 



each of the /%(s) 



1,2, j = 0,1,2, into the sum of an overall term and 



a local adjustment. That is, 



(5) 



The six-dimensional process (/3ij(s))j = i 5 2 ; j=o,i,2 is in turn modeled cLS 3> SIX - 
dimensional mean-zero correlated Gaussian process, again using the method 
of coregionalization. Therefore, we express each of the spatially-varying co- 
efficients Pij(s) as a linear combination of mean-zero unit-variance latent 
independent Gaussian processes, each equipped with an exponential covari- 
ance structure. More specifically, assuming, without loss of generality, that 
the coregionalization matrix A is lower triangular, we have that: 



(6) 



/ fto(s) \ 
fti (s) 

012 (S) 
020 (S) 
021 (S) 
V 022 (S) J 



( wi(s) \ 
w 2 (s) 
w 3 (s) 

104 (S) 

10 5 (S) 



where, for each k = 1, . . . , 6, Cov(iOfc(s), Wk(s')) = exp(- 



9fe s- 



with <j>k 



decay parameter for i0fc(s). As in the univariate case, the coregionalization 
matrix A determines, through T = AA', the correlation between the six 
spatially- varying coefficients 0ij(s). Additionally, as a consequence of (4), A 
induces a correlation structure on the bivariate random vector Y(s), i.e. 



(7) E 



Y(s),Y(s') 



[I 2 ® (1 x\{B) x 2 (B))'] AS^A'- [1 2 ® (1 x x {B') x 2 {B'))} 



where B and B' are grid cells containing, respectively, s and s', I2 is the 
identity matrix of order 2, and S m is a 6 x 6 diagonal matrix with i-th 
element exp(— ^||s — s'||), i = 1, . . . , 6. 

Our bivariate downscaler model, though simple, is general and flexible. 
The specification of the bias of the numerical model by means of spatially- 
varying coefficients recognizes that calibration under the numerical model 
should be site-specific. Again, as in the univariate case, a constant-mean bi- 
variate downscaler can be obtained as a special case. Moreover, permutation 
of the entries in the 0(s) vector does not affect the prior. The joint model 
still presents a 6 x 6 local covariance matrix modeled through its Cholesky 
decomposition. 

Simplified versions of the general bivariate downscaler can be obtained by 
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setting entries in the coregionalization matrix A equal to zero. 2 For example, 
if the set of non-null entries in the A matrix is given by {Ai, A21, A22, Am, Aj4, 
and we assume that /3i2 = and fat = 0, then our bivariate downscaler 
reduces to two independent univariate downscalers on ozone and PM2.5, re- 
spectively. This form of A implies that the spatially- varying coefficients are 
correlated within pollutants but not across pollutants. 

A simpler version of the bivariate downscaler that requires only three 
non-null entries in the coregionalization matrix A but still induces a cor- 
relation among the two components Vi(s) and Y2(s) of Y(s) is the model 
that assumes that only the local intercept adjustments, /?io(s) and /?2o(s), 
are non-null and correlated. In this case, the coregionalization matrix A has 
only three non-null entries, {An, A41, A44} and the covariance between the 
square root of the observed ozone concentration at s and the logarithm of 
the observed PM2.5 concentration at s' reduces simply to: 

(8) Cov^s),!^')) = A 1 i-^4iexp(-0 1 ||s-s , ||) 

A simple extension to this model can be achieved by maintaining the 
same correlation structure for the bivariate random vector Y(s), that is, by 
assuming again that An is the only non-null off-diagonal element of the core- 
gionalization matrix A, but by postulating that the local adjustments to the 
coefficients of the two numerical model outputs, /3n(s), /3i2(s), /32i(s), and 
/?22(s), are non-null independent Gaussian processes with variances equal, 
respectively, to A^, A33, A\^ and A\§. Then, the coregionalization matrix 
A corresponding to this model has seven non-null entries: An, A22, A3, 
An, Am, A55, and A m . 

An extension of this last model that still involves fewer parameters than 
the general formulation for the bivariate downscaler (13 versus 21) is the 
model which assumes that only local adjustments of the form j3ij(s) and 
fl2j(s), for j = 0, 1,2, are correlated across pollutants, while all the /3y(s)'s 
specific to the same pollutant are correlated. Such a model specification 
corresponds to the following coregionalization matrix: 



(9) A = 



/ A11 














\ 


An 


A22 














A31 





A33 











Ai 








Am 











A52 





A34 


A35 





V 





A33 


A34 





Am J 



2 In principle, we could include all 21 Aij's in the model and let the data suggest which 
are significant. Instead, we have chosen to do model comparison between several models, 
each having plausible interpretations. 
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Other simplifications of the bivariate downscaler general model can be en- 
visioned. 

The specification of our Bayesian hierarchical bivariate downscaler model 
is completed with priors on the parameters. We adopt an inverse gamma 
with large variance for the nugget variances, rf and r|; independent vague 
normals for the overall regression coefficients f3ij, i = 1,2, j = 0,1,2; log- 
normals with vague standard deviations for the diagonal entries of the core- 
gionalization matrix A; and vague normals for the off-diagonal elements of 
A. Details on handling the decay parameters k = 1, . . . , 6, are provided 
in Section 4.1. 

3.3. Space-time bivariate downscaler. We now extend our general bivari- 
ate downscaler to accommodate data that has been collected over time. Let 
t = 1 , . . . , T denote the times at which we have observations and numerical 
model outputs, and let Yi(s, t) and Y2(s,i) denote, respectively, the square 
root and the logarithm of the observed ozone and PM2.5 concentration at 
site s on day t. In an analogous way, let x\{B,t) and X2(B,t) indicate, re- 
spectively, the square root and the logarithm of the CMAQ output for ozone 
and PM2.5 on day t over the 12-km grid cell B. Then, if s lies in grid cell B, 
the model for our bivariate downscaler becomes 

Yi(s,t) = f3 1 o(s,t)+Pi 1 (s,t)xi{B,t)+j3 1 2(s,t)x 2 (B,t) + ei{s,t) 

(10) y 2 (s,t) = p 2 o(s,t)+f32i(s,t)x 1 (B,t) + i322(s,t)x2(B,t) + e 2 (s,t) 

where ei(s,i) and e 2 (s,t) are two white noise processes that follow indepen- 
dently a N(0, t\ ) and a N(0, r| ) distribution. 

As in the spatial setting, for each i = 1, 2 and j = 0, 1, 2, we write: 

(11) hjM = Pij,t + PijM 

where, for each t = 1, . . . ,T, [3ij (s, t) are correlated Gaussian processes. 

Following [3], we can model the temporal dependence in the components 
of (11) in several ways. For example, for the f3ij/s, we can assume that they 

are either independent across time, i.e. /3ij s t ^ iV^/i^, cr?), or, alternatively, 
that they evolve dynamically in time [52], that is, 

(12) p ijJt = (Hjfat-i + mj,t, Vij,t N(0, 4), % = 1, 2; j = 0, 1, 2 

and /% j0 ~ ^(ftj,o,4o)' 

Similar time-dependence can be imposed on the Gaussian processes (3ij(s, t). 
We can assume that the correlated Gaussian processes f3ij(s,t) are of the 
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form 



(13) 



/ 0ioM) \ 
0ii (M) 
/MM) 

020 (M) 

021 (S,t) 

V 022 (S,t) / 



A 



/ u>i(s,i) \ 

102 (S,t) 

10 3 (M) 
Wi(s,t) 

w 5 (s,t) 



V w 6 (s,t) / 



with A still a coregionalization matrix and with the underlying Gaussian 
processes, Wk(s,t), nested within time. Then, at each t = 1,...,T, the 
u>fc(s, t) are independent replicates of mean-zero unit- variance Gaussian pro- 
cesses with exponential correlation and decay parameters 4>k,t that can be 
either taken to be constant in time or independent across time. 

In the second case, we can model the local adjustments, 0ij(s, t), to evolve 
dynamically in time. As in [17], for each i = 1, 2 and j = 0, 1, 2, we assume 
that 



where the innovations Vij(s,t) are correlated Gaussian processes. In other 
words, in this second case, the coregionalization (13) is specified on the 
Vij(s, t), rather than on the 0y(s, t), but it still employs the underlying Gaus- 
sian processes Wk(s,t), which are defined as above. This dynamic model for 
the 0ij(s, t) is completed by specifying as initial conditions that 0jj(s, 0) = 
for z = l,2 and j = 0, 1, 2. In both cases we can be more general and assume 
that the coregionalization matrix in (13) is indexed by time, with entries 
independent across time. Furthermore, if t is on a continuous domain, then 
the Wk(s,t) can be independent space time processes. 

The two different time dependence structures for the overall regression 
coefficients, 0jj,t, and for their local adjustments, 0ij(s, t), can be combined 
together in four ways, yielding four models that correspond to different as- 
sumptions on the way the overall and local performance and calibration of 
the numerical model changes over time. 

3.4. The general multivariate downscaler. Extension to a multivariate 
downscaler given observational data and numerical model outputs for p ran- 
dom variables is evident. As in the bivariate case, we start with the static 
formulation of the model and then we extend it to the spatio-temporal set- 
ting. 

Let Yi(s), i = 1, . . . ,p be the observed data for the i-th variable at a site s 
in the spatial domain S and let Xi(B), i = 1, . . . ,p, be the numerical model 



(14) 
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output for the i-th variable over grid cell B. Then, again we associate to 
each site s the grid cell B in which s lies and we relate the observational 
data and the numerical model output as follows: 

(15) Yi(s) = Ao(s) + f]4-(s)x i ( J B) + e i (s) e^s)^ N(0,t?) 

i=i 

for i = 1, . . . ,p. 

We decompose each of the p ■ (p + 1) terms /3{j(s), i = 1, . . . ,p, j = 
0,1, ... ,p in the sum of an overall term and a local adjustment: Pij(s) = 
Pij + ftj(s), and we model the /3j,-(s)'s as correlated mean-zero Gaussian 
processes with exponential covariance structures using again the method of 
coregionalization. In the general p-dimensional multivariate downscaler, the 
coregionalization matrix A is a p(p + 1) x p(p + 1) matrix. Specifications 
of (15) similar to the above can be considered by simply setting to zero 
entries of the A matrix, thus inducing simplifications on both the correlation 
structure of the multivariate random vector Y = {(ii(s))i = i p : s 6 S} and 
on the covariance structure of the single Yi(s), i = 1 . . . ,p. 

To complete the specification of the model, we place standard priors on the 
model parameters: thus, we specify p(p + 1) vague normals for the overall 
regression coefficients fy; p inverse gammas with large variances for the 
nugget variances rf, i = 1 . . . ,p; log-normals with large variances for the 
diagonal entries A^k, k = 1, . . . ,p(p+ 1), of the coregionalization matrix A; 
and vague normals for the off-diagonal entries of A. 

If the data on the p random variables has been collected not only over 
space, but also over time, then the multivariate general downscaler would 
be extended to a spatio-temporal setting in an obvious way. Let Yi(s,t), 
i = 1, ... ,p be the observed data for the i-th random variable at site s at 
time t, t = 1, . . . , T, and let Xi(B, t) be the model output for the i-th variable 
over grid cell B at time t. Then, for i = 1, . . . ,p, we postulate the following 
relationship between observational data and numerical model output: 

p ~ -a 

(16) Yi(s,t) =p i0 (s,t)+J2Pij(s,t)x j (B,t)+e i (s,t) e;(s,t) ~ N(0,rf) 

3=1 

where s lies in grid cell B. 

For each i = l,...,p and j = 0,1,..., p, we write: (3ij(s,t) = fiijj + 
(3ij(s, t), and we model the temporal structure in the and in the /3jj(s, t), 
following what we have presented in Section 3.3 for the bivariate downscaler. 
Thus, we assume that the fiij/s are either nested within time or dynamic in 
time, and, similarly, we postulate that the fiij{s,t) are either independent 
replicates over time or are dynamic in time. 
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4. Model fitting. 

4.1. Priors. We have already mentioned briefly in Sections 3.2 and 3.4 
the form of the prior distributions for the parameters of the bivariate and 
multivariate downscaler both in their static and spatio-temporal formula- 
tions. Here, we turn to the spatial decay parameters (ftk of the latent Gaussian 
processes wj~(s) used in the coregionalization. Discrete uniform priors facili- 
tate model fitting. Previous experience with such priors has always resulted 
in posterior distributions placing highest probability on a value which is, 
essentially, the Restricted Maximum Likelihood (REML; [24, 38]) estimate. 
Hence, in what follows, we propose to keep the spatial decay parameters 
fixed and equal to values that are determined by a sensitivity analysis. 

More specifically, consider the simplified bivariate downscaler that uses 
only the two latent Gaussian processes wi(s) and w^(s) and is obtained when 
the coregionalization matrix A has only three non-null entries, An,A4i, 
and A44. The decay parameters <pi and $4 determine the way the spatial 
correlation decays with distance in w\{s) and W4 (s), respectively, and, by 
consequence, in Yi(s) and ^(s). In fact, the covariance structure of i^(s) 
induced by such simplified version of the bivariate downscaler is given by 
the sum of two exponential covariance functions with decay parameters, 
respectively, equal to and $4, and a diagonal matrix with entries all 
equal to the nugget variance r| • 

To obtain a rough estimate of the magnitude of the decay parameters 
01 and (j>4, we have proceeded as follows. We have considered the moni- 
toring data for the square root of ozone and the logarithm of PM2.5 for a 
given day t. We have modeled the monitoring data for both variables as 
Gaussian processes with an unknown mean, respectively, fj,y lt and /Uy 2 t , 
and with an exponential covariance function plus nugget effect, with decay 
parameters 0y x t and </>y 2t , and nugget variance, t y t and Ty 2 ■ We have esti- 
mated the parameters via Restricted Maximum Likelihood (REML; [24, 38] 
as implemented in the geoR package of R. We have repeated this operation 
for each day in the period June 1- September 30, 2002 and have obtained 
daily estimates of 4>y 1 t and 4>Y 2 1 ■ Since histograms of those daily estimates 
exhibited long right tails, we have summarized those distributions by tak- 
ing the median of the daily estimates, which yielded, respectively, 0.0016 
and 0.00125, corresponding to ranges of, respectively, 1875 and 2400 km. 
We have then performed a sensitivity analysis to determine how the pre- 
dictive performance of a downscaler model is affected by the magnitude of 
the decay parameter. Therefore, keeping the decay parameters 4>l anci 4>i 
fixed and equal to the estimated medians of §y x t and 0y 2 1 , we have fit the 
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spatio-temporal version of the bivariate downscaler that models the overall 
terms /%,t, i = 1,2, j = 0,1,2 and the local adjustments to the overall 
intercepts, /3jo(s,t), i = 1,2 as independent in time (see Section 7.2 in the 
Supplemental material online). For each day, we have predicted levels of 
ozone and PM2.5 concentration at the validation sites so by sampling from 
the posterior predictive distributions, /(Y(so)| {Y(s)} , {x(5)}) and subse- 
quently backtransforming the prediction to the original scale. We have then 
compared the predictions with the observations by computing the predictive 
mean square error (PMSE), predictive mean absolute error (PMAE), empir- 
ical coverage of the 95% predictive interval and width of the 95% predictive 
interval. More specifically, denoting with Y*(s,t) the backtransformed data 
(recall we used the square root and the logarithm transform on ozone and 
PM2.5, respectively), we have computed the PMSE as follows: 

T V 

(17) PMSE = — (Y*{ar,t) - Y*(s r ,t)) 2 • I(F*(s r ,i)) 

nr " t=\ r=l 

where n v denotes the total number of observations in the validation dataset 
(6530 for ozone or 2559 for PM2.5 as noted in Section 2), s r is the r-th site in 
the validation set containing a total of V sites, Y*(s r , t) denotes the posterior 
mean of the predictive distribution at site s r at time t on the original scale, 
and I(Y*(s r ,i)) is equal to 1 if Y*(s r ,t) is observed and otherwise. An 
analogous definition holds for the PMAE, with Y*(s r ,t) now referring to 
the posterior median of the predictive distributions, again on the original 
scale. We have generated predictions and computed the summary statistics 
mentioned above four times, each time keeping the decay parameters fixed 
and setting them equal, respectively, to the estimated medians of 1 and 
0y 2 t , multiplied by 10 and 100, and divided by 10 and 100. 

To compare the performance of the bivariate downscaler to that of a uni- 
variate downscaler, we have performed the same sensitivity analysis also for 
the univariate downscaler. Thus, keeping the decay parameters for ozone 
and PM2.5 fixed and equal to the values considered for the bivariate down- 
scaler, we have fit two spatio-temporal univariate independent downscalers 
for ozone and particulate matter, where the only non-null local adjustment 
term was the adjustment to the overall intercept, and where both the over- 
all regression terms and the spatially varying coefficients were nested within 
time. Then, as in the bivariate downscaler case, we have predicted ozone and 
PM2.5 concentration at the validation sites and we have assessed the per- 
formance of the predictions using the same summary statistics mentioned 
above. 

We have carried out a similar sensitivity analysis also for ordinary kriging 
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[8, 9]. Therefore, for both ozone and PM2.5, for each day, using the same 
range of values for the decay parameters used in the experiments with the 
univariate and bivariate downscalers, we have kriged the observed data, on 
the transformed scale, to the validation sites, and we have subsequently back- 
transformed the predictions. Finally, we have compared the predictions with 
the observations using the same summary statistics that we have employed 
for both downscalers. 

Table 1 and Table 2 report, respectively, the Predictive Mean Square Er- 
ror (PMSE), Predictive Mean Absolute Error (PMAE), the empirical cov- 
erage and the width of the nominal 95% predictive interval for ozone and 
PM2.5 for the two space-time downscalers and for ordinary kriging. In both 
tables, the results corresponding to predictions obtained when the models 
have been fitted using values of the decay parameters equal to the medians 
of the daily estimates of 4>y 1 t and 4>y 2 1 are reported in the middle column. 
From both tables, it is clear that when the decay parameter is one and, 
especially, two orders of magnitude larger, the quality of the predictions ob- 
tained using both downscalers deteriorates noticeably in terms of PMSE and 
PMAE. Also, the predictive intervals are wider in this case, an indication 
that there is increased variability in the predictions. Misspecifying the mag- 
nitude of the decay parameters on the small side affects the quality of the 
predictions less noticeably. So, altogether, the sensitivity analysis has shown 
that both the univariate and the bivariate downscaler and ordinary kriging 
yield predictions that have overall best validation statistics when the decay 
parameters are equal to the median of the daily estimates of 4>y 1 1 and <f>Y 2 1 
and so, in the sequel, we fix them at these values. 

4.2. Handling misalignment. As noted in Section 2, not all sites report- 
ing ozone concentration measure PM2.5 and vice versa. This creates a spatial 
misalignment in the data that we handle through the latent Gaussian pro- 
cesses Wk(s) involved in the coregionalization, after having appropriately 
permuted and partitioned the data vector. Here, we illustrate how to han- 
dle the spatial misalignment in the static setting. Extension to the spatio- 
temporal case is straightforward. 

Let St be the set of nt sites reporting measurements of ozone and/or 
PM2.5 concentration on day t. We can decompose St as the union of three 
disjoint sets, 5 Bo th,t, S ,t, and S PM ,t- 

(18) S t = 5Both,t U So,t U SpM,t 

where the first set includes all sites s £ S in which both ozone and PM2.5 
have been measured on day t, the second contains all the sites in which only 
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ozone was measured on day t, and the last includes all sites s where only 
PM2.5 has been measured on day t. 

Following the decomposition in (18), we reorder and partition the random 
vectors Y = (Yi(s), l2(s))seS< and = (wfc(s)) se s t in the following way: 
Y = (Y Bo th, Y , Y PM ) and w fc = (w fe)Both , w fc>0 , w k)PM ) where Y Bo th = 
(Yi(s), y 2 (s))se5 Both t and analogous definitions hold for Y , Y PM , w fc>Both , 
w ki0 , and w k}PM . 

Then it is clear that both components Yi ]Bo th = (Yi(s)) s& s B , and 
Y2, Bo th = (Y2(s)) sg 5 B h t have non-missing values and both contribute to 
the log-likelihood (??), while for s 6 So,t and s E SpM,t only one component, 
respectively, Y^o and Y2 i pm, has non- missing values and contributes to the 
log-likelihood. 

The latent Gaussian processes w k (s) are defined on the entire spatial do- 
main S. To obtain samples from the posterior distribution of w k (&), within 
each MCMC iteration, we draw samples from the full conditionals by pro- 
ceeding in the following way. If G denotes the collection of all model pa- 
rameters (thus, G = (Ti,T$,{A k l} k>l=li ...,6;i<fc!{/%}i=l,2;j=0,l,2))' for each 

k = 1, . . . , 6, we: 

i. sample w fc)Both from the full conditional 7r(w fcj Both|w fct o, ™k,PM, Y Bot h, G); 

ii. sample w kj o from the full conditional vr(wfc 5 o| w fc,Both, Wk,PM, Y^q, Q); 

iii. sample w kjP M from the full conditional ^(w^pmI w^Both, w fc,0) Y2,?m, ©); 

iv. for each s E S \ St, sample w k (s) from the conditional distribution 

n(w k (s)\Wk ) Both,VV r k) 0,Wk > PM, ©) 

By repeating steps i-iv within each MCMC iteration, we obtain samples 
from the posterior distribution of w k (s) for s E S, thus solving the problem 
of spatial misalignment between ozone and PM2.5. 

4.3. Model selection. Using results from an analysis on the predictive 
performance of the univariate spatio-temporal downscaler for ozone concen- 
tration [3], we consider only the spatio-temporal versions of the bivariate 
downscaler with time-varying parameters independent across time. More 
specifically, we have fitted the spatio-temporal versions of the four bivariate 
downscalers that we have presented in Section 3.2, that is: 

i. the bivariate downscaler equivalent to two independent univariate down- 
scalers, obtained when the coregionalization matrix A has as non-null 
entries only 

{An, A21, A22, A44, Aq4, Am} ', 
ii. the bivariate downscaler with only {.An, A41, A44} as non-null entries 
in the coregionalization matrix A; 



imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: April 8, 2010 



A BIVARIATE SPACE- TIME DOWNSCALER 



19 



iii. the bivariate downscaler with {An, A22, A33, A41, A44, A55, ^66} as non- 
null entries in the coregionalization matrix A; and 

iv. the bivariate downscaler with coregionalization matrix A given by (9) 

For each of these models, we have examined their out-of-sample predic- 
tive performance, by predicting daily concentrations of ozone and PM2.5 
at the 65 validation sites described in Section 2 and shown in Figure 1. 
These predictions are compared with the observed ozone and PM2.5 lev- 
els at the validation sites during the period June 1-September 30, 2002. 
Note that, for each day, predictions of ozone and PM2.5 at a validation 
site so are obtained by sampling from the posterior predictive distribu- 
tion /(Y(so, t)\ {Y(s, t)} , {x(i?, t)}) and then by subsequently transforming 
them back to the original scale. 

To quantify the predictive performance of each bivariate downscaler we 
have utilized the same validation statistics used in our sensitivity analysis 
presented in Section 4.1. We also employ two proper scoring rules: (i) the 
Continuous Ranked Probability Score (CRPS; [21]) defined as: 



where F is the cumulative predictive distribution function, y is the observa- 
tion that materializes, 1 is the Heaviside function, that is equal to 1 if z is 
greater than y and otherwise, and (ii) the Interval Score [21], defined for 
a (1 — a) ■ 100% predictive interval with lower bound I and upper bound u 
as: 



where y denotes again the observed value. 

So, for each model, we have computed the Predictive Mean Square Error 
(PMSE), Predictive Mean Absolute Error (PMAE), the empirical coverage 
of the 95% predictive interval, and the width of the 95% predictive inter- 
val to obtain an indication not only of the average bias and errors in the 
predictions, but also of the level of uncertainty in the predictions. Addition- 
ally, the CRPS, being a strictly proper scoring rule, provides a simultaneous 
assessment of the calibration and sharpness of the posterior predictive dis- 
tribution whereas the Interval Score rewards predictive distributions with 
narrower predictive intervals while imposing a penalty for observations lying 
outside the predictive interval. 

Finally, to determine the improvement in the predictions of ozone and 
PM2.5 obtained from using our downscaling approach with respect to non- 
model based techniques, we have compared the downscaler with ordinary 
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kriging and cokriging [8, 9]. For kriging, using as decay parameters for ozone 
and PM2.5 concentration, respectively, 0.0016 and 0.00125, for each day, 
we have kriged, separately, the observed daily concentrations of ozone and 
PM2.5 at the test sites to the validation sites, working on the transformed 
scale and then backtransforming the predictions. For cokriging, we have ex- 
ploited the information contained in the copollutant and, for each day, we 
have cokriged the observed daily concentrations of ozone, first, and PM2.5, 
afterwards, to the validation sites. Note that in performing cokriging, we 
had to define, not only a covariance function for both ozone and PM2.5, but 
also a cross-covariance function. For both pollutants, as covariance function 
we used again the exponential covariance function with decay parameters, 
respectively, equal to 0.0016 and 0.00125. For cross-covariance function, in- 
stead, we used the cross-covariance function induced by the second version of 
the bivariate downscaler listed above, explicitly given in (8). The parameters 
for such cross-covariance function were set to be equal to the posterior means 
of j4ii,v44i and A44 obtained from fitting the bivariate downscaler model in 
ii. Finally, predictive intervals and uncertainty for kriging and cokriging are 
based on usual formulas for the kriging and cokriging variance [8]. 

5. Data Analysis. We first present validation results for the different 
models. Since we did not find a significant improvement in the predictions 
of ozone and PM2.5 obtained from versions iii and iv of the bivariate down- 
scaler, we show only results for the versions, i and ii. To distinguish between 
the two, we will refer to the first as the independent downscaler, while the 
second will be called the bivariate downscaler. Comparing the performance 
of the two models allows us to establish the gain in predictive performance 
that can be ascribed to explicitly accounting for the correlation between the 
two pollutants. 

From Section 4.3, for each of the models considered, we have computed 
the predictive mean square, the predictive mean absolute error, the empirical 
coverage, the width of the 95% predictive intervals, the Continuous Ranked 
Probability Score and the Interval Score, averaging across sites and days. In 
order to determine whether distance from the closest monitoring sites has 
an effect on the quality of the predictions, we have divided the validation 
sites into two groups: those that have a distance from the closest monitoring 
test site of less or equal than 40 kilometers, and those that are more than 40 
kilometers far apart from the closest monitoring test sites. For each of the 
two groups, consisting, respectively, of 13 and 52 sites, we have computed 
the same validation statistics mentioned above. 

Table 3, Table 4, and Table 5 report all these summary statistics, respec- 
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tively, for all the validation sites, for the validation sites that are less or 
equal than 40 km distant from the closest monitoring test sites, and for the 
validation sites that are more than 40 km away from the closest monitoring 
test site. 

From all tables, it is clear that the bivariate downscaler yields predictions 
of both ozone and PM2.5 concentrations that are less biased than those 
obtained using the univariate downscaler. Additionally, the bivariate down- 
scaler has a lower average CRPS score, indicating that its predictive distri- 
bution is sharper and better calibrated than the one corresponding to the 
univariate downscaler. In terms of width and empirical coverage of the 95% 
predictive intervals, both downscalers perform similarly and both have em- 
pirical coverage close to nominal. However, the combined assesment of both 
properties, which is provided by the interval score, favors again the bivari- 
ate downscaler over the univariate downscaler. A similar trend in predictive 
performance can be observed for kriging and cokriging: exploiting the corre- 
lation between ozone and particulate matter yield predictions that are less 
biased than those obtained by not accounting for it. Additionally, kriging 
always underestimates uncertainty in the predictions, producing predictive 
intervals that are too narrow and, thus, do not have the nominal coverage. 
This is also reflected in the CRPS and the Interval Scores: kriging always has 
a lower CRPS score than cokriging probably due to the smaller widths of the 
predictive intervals, which in turn, failing to achieve the nominal coverage 
are penalized in terms of Interval Score. 

Comparing the bivariate downscaler with cokriging, we can see that the 
two methods perform equally well in terms of Interval Score; however the 
former outperforms the latter in terms of PMSE, PMAE, and CRPS for both 
ozone and PM2.5. This follows from the fact that the bivariate downscaler 
not only models the correlation between ozone and PM2.5 but also uses the 
information contained in the numerical model output, while cokriging does 
not take advantage of this additional information when predicting levels of 
ozone and PM at the validation sites. In addition, the bivariate downscaler 
predicts levels of PM2.5 better than cokriging, likely because of the sampling 
frequency of PM2.5 which renders it rather difficult to interpolate to the en- 
tire spatial domain, especially when the number of sites with observations 
is fairly low. Note that validation of the out-of-sample predictive perfor- 
mance of both downscalers and kriging/ cokriging for PM2.5 occurs mostly 
every three days due to the sampling frequency of PM2.5. In days with fewer 
observations from monitoring sites, this can diminish predictive gain associ- 
ated with both downscalers relative to spatial interpolation techniques based 
solely on monitoring data. 
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The separation of validation sites into two groups based on their distance 
to the closest monitoring test site, does not reveal a difference in terms of the 
predictive performance of the various models. The bivariate downscaler still 
outperforms all the other models. However, when considering sites that are 
less than 40 kilometers far from the closest monitoring test site, the bivariate 
downscaler and cokriging predict ozone concentration equally well. This is 
not true for particulate matter, nor it is true when we consider validation 
sites that are more than 40 kilometers away from the closest monitoring 
test site. For this group of sites, for both ozone and PM2.5, the predic- 
tions obtained using the bivariate downscaler are superior to the predictions 
obtained using any other method. Counterintuitively, for all methods, the 
quality of the predictions are better at validation sites that are farther from 
the closest monitoring test sites than at validation sites that are closer to 
the monitoring test sites. 

Our bivariate downscaler model, as well as the univariate version, can 
also be used to generate predictive surfaces for both ozone and PM2.5 con- 
centration. Figure 5 and Figure 6 display, in panels (d), predictive surfaces 
of ozone and PM2.5 concentration for June 25, 2002 obtained using our bi- 
variate downscaler model. Each panel displays the posterior predictive mean 
of ozone and PM2.5 concentration obtained by sampling from the posterior 
predictive distribution at each site so on the CMAQ grid covering the study 
region. As a comparison, Figure 5 and Figure 6 present also, in panel (c), 
the predictive surfaces obtained by simply kriging the observed concentra- 
tions of ozone and PM2.5 to the 10,504 CMAQ grid cells covering the region. 
Finally, panels (a) and (b) show respectively the observed and the CMAQ 
predicted concentrations for both pollutants. 

As both figures show, the predictive surfaces obtained using kriging are 
very smooth for both pollutants, while the surfaces obtained using our bi- 
variate downscaler show more texture which is expected to provide improved 
prediction. In particular, our bivariate downscaler method generates predic- 
tive surfaces that correct the CMAQ predictions for local bias, and display 
a gradient that reproduces some of the features of the CMAQ surfaces while 
accounting for the spatial gradient in the observations. For example, Fig- 
ure 5(d) predicts lower level of ozone than CMAQ in the South-Eastern 
corner of map, while Figure 6(d) pushes up the level of PM2.5 predicted by 
CMAQ in the same region. 

In order to visually quantify the local biases of the CMAQ predicted 
concentrations of ozone and PM2.5 for June 25, 2002, Figure 7 displays 
spatial maps of the posterior predictive means of Pio(s) and /?2o( s )- Since 
the bivariate downscaler models has been developed on the transformed 
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scale, immediate conversion of these maps into calibration for the numerical 
model is not possible. However, they do display the spatial variability in 
the additive bias of the numerical model output as well as identify areas 
where the numerical model either underpredicts or overpredicts a pollutant 
concentration. 

An additional feature of our downscalers is the possibility to upscale and 
predict the variable of interest at any spatial scale. To illustrate such capabil- 
ity for our bivariate downscaler, we have divided our study region into three 
subregions: the Mid-Atlantic (Virginia, West Virginia, Maryland, Delaware, 
and Pennsylvania), the Midwest (Kentucky, Ohio, Illinois, Indiana, and Mis- 
souri), and the South (all the remaining states contained in the study re- 
gion). For each of these regions, we have predicted the average ozone and 
PM2.5 concentration on June 25, 2002, by predicting on a grid of points us- 
ing the posterior predictive distribution induced by our bivariate downscaler 
model and by averaging the predictions over each of the three regions. Those 
predictions were then used to consider contrasts that might be of interest. 
Hence, Figure 8 presents plots of the posterior predictive distribution of the 
contrast in average ozone and PM2.5 concentration, respectively, between the 
Mid- Atlantic and the Midwest, Mid- Atlantic and the South and the Midwest 
and the South for June 25, 2002. As Figure 8 shows, both the average ozone 
and PM2.5 concentration are higher in the Mid-Atlantic region than in the 
Midwest and the South on June 25, 2002, with the South having, on average, 
the lowest average concentration for both pollutants on the selected day. We 
could further push such comparisons to look at averages or contrasts over 
time. 

6. Summary. We have proposed a bivariate downscaler model that al- 
lows us to scale down the outputs of a numerical model from grid level 
to point level. The model is rather general and flexible and can easily be 
extended to a multivariate downscaler model for p numerical model out- 
puts. For the bivariate case, we have examined several simplifications of the 
general bivariate downscaler model, each inducing a different type of corre- 
lation structure among the two components of the bivariate random vector. 
We have also shown how the model can be extended to accommodate data 
collected over both space and time. Also, as a process model, by suitable 
block averaging, inference can be scaled up to provide average exposure over 
arbitrary regions. In addition, we have demonstrated the computational ad- 
vantages of scaling down, due to the relative sparsity of monitoring data 
compared with model output data, and the computational feasibility for 
our models with large space-time datasets. We acknowledge that the small 
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decay parameters (equivalently, large range parameters) relative to the size 
of our region, suggest that a linear variogram associated with a valid co- 
variance function over a bounded region might be an alternative modeling 
path. In turn, this suggests the possibility of an approach based entirely on 
contrasts, though, given that our interest is in spatially varying regressions 
on the model output, it remains unclear how to proceed if such a modeling 
choice was undertaken. 

In our case study, we have applied the bivariate downscaler to downscale 
and predict ozone and PM2.5 concentration during the summer months of 
year 2002; however it is clear that our model can be used for any pair 
of pollutants. For example, it would be interesting to apply the bivariate 
downscaler model to predict wet and dry sulfate deposition, embedding the 
bivariate downscaling approach in the modeling context proposed by [42]. 
Another potentially interesting analysis would be to work with O3 and CO. 

There are two extensions we are currently exploring. First, in both the 
univariate and the bivariate downscaler model, the numerical model output 
has been taken as data and we have not accounted for uncertainty in the 
predictions given by the numerical model. Here, potential issues involve too 
much smoothing of the CMAQ output, potential displacement of CMAQ 
grid cells, and errors in inputs that drive CMAQ. Within our downscaling 
approach it is possible to characterize the uncertainty in the numerical model 
output by adding another hierarchical level to our framework. That is, we 
return to the issue in the Introduction where we now require a stochastic 
model for the CMAQ output. 

Secondly, we seek extension beyond Gaussian specifications. For instance, 
if we were to study extreme fields for two pollutants from both monitor- 
ing data and CMAQ, we would not use Gaussian models. Similarly, if we 
recorded two binary variables at each location, resulting in a 2 x 2 table for 
the location, again, Gaussian processes are not appropriate. 

7. Appendix. 

7.1. Covariance and cross- covariance. Here we present explicit formulas 
for the covariance and cross-covariance between Yi(s) and Y2( s ) for the 
different versions of the bivariate downscaler considered in Section 3.2. 

For the bivariate downscaler with coregionalization matrix A with non- 
null entries, An, A41, A44, the covariance between the square root of the 
observed ozone concentration, Y\ (s) and the logarithm of the observed PM2.5 
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concentration at s' is given by (8). Furthermore, for this model: 



Cov(yi(s),yi(s')) 
Cov(y 2 (s),y 2 (s / )) 



A? 1 exp(-^>i||s- s'||) +rf J SiS / 

A 2 Al exp(-^i||s - s'||) + ^4 4 exp(-0 4 ||s - s'||) + r|(5 s 



where <5 SjS / is equal to 1 if s = s' and equal to otherwise. 

For the bivariate downscaler with coregionalization matrix A with non- 
null entries, {An, A22, -A33, -A41, .A44, ^55, ^66}, the cross-covariance between 
Y\ (s) and ^(s') is given by (8), while the expressions for the within pollutant 
covariances are, respectively: 



Cav(Yi(s),yi(s / )) 



An exp(— 0i||s — s'| 
3 



+ [ A lk ■ Xfc_i( J B)x fe _i( J B / ) • exp(-<^||s - s'||)] + t?6 b 



k=2 



Cov(F 2 (s),y 2 (s')) = Al 1 exp(-</) 1 ||s-s , ||) +^ 4 exp(-0 4 ||s-s'| 



^ [ A kk ■ x {k-4){ B ) x {k-4){ B ') ■ exp(-0 fc ||s - s'| 
fc=5 



2 u s,s' 



with B and B' grid cell containing, respectively, s and s'. 

Finally, for the bivariate downscaler model with coregionalization matrix 
A given by (9), the covariance between the two components Y\(s) and ^(s') 
of the bivariate random vector Y(s) is given by: 



Cov(y 1 (s),y 2 (s / )) 

(19) 



= AxxMi • exp(— 0i||s — s'||) 
3 

+ Y [ A kkA( k+3 ) k • x k -i(B)x k -x(B') • exp(-^ fc ||s - s'| 

k=2 



while the expression for the inter-pollutant covariances are given, respec- 
tively, by: 

Cov(Yx{s),Yx(s')) = AneM-Ms-s , \\)-[An + A2 1 (xx(B)+x 1 (B')) 

+ A 3 x(x 2 (B) + x 2 (B'))] 
3 

+ Yl [ A lk ■ x k _ 1 (B)x k ^ 1 {B') ■ exp(-<^||s - s'||)] + rfo 

k=2 
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and 



Cov(Y 2 (s),r 2 (s')) 



A\ x exp(-^>i||s - s'||) 



3 



+ 



S [ A lk+S)k x k-i{B)x k - 1 (B') • exp(-0 fc ||s - s'||) 



fc=2 



+ 



+ 



A 44 • exp(-0 4 ||s - s'||) • [A 44 + A 54 ( Xl (B) + xi(B')) 
A 64 (x 2 (5) +x 2 (#'))] 



6 



+ 



A lk • Zfc-4(-B)zfc-4(-B') • exp(-^ fe ||s - s' ||) + t|<5 SjS / 



fc=5 



Disclaimer: 

The U.S. Environmental Protection Agency through its Office of Research 
and Development partially collaborated in the research described here. Al- 
though it has been reviewed by the Agency and approved for publication, it 
does not necessarily reflect the Agency's policies or views. 
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Table 1 

Predictive Mean Square Error, Predictive Mean Absolute Error, empirical coverage and 
width of the 95% predictive interval for O3 for different values of the decay parameter 

under different models. 



Model 


Summary 
statistic 


Decay parameter 


0.16e-4 


0.16e-3 


0.16e-2 


0.16e-l 


0.16 


Independent 
downscaler 


PMSE 


55.54 


55.94 


55.92 


58.13 


95.17 


PMAE 


5.54 


5.56 


5.56 


5.68 


7.39 


Emp. cov. 


93.2% 


92.0% 


92.9% 


91.7% 


92.9% 


Width 


27.33 


26.30 


27.34 


27.19 


35.78 


Bivariate 
downscaler 


PMSE 


53.21 


51.64 


52.78 


55.24 


91.42 


PMAE 


5.45 


5.36 


5.40 


5.54 


7.27 


Emp. cov. 


92.6% 


92.3% 


93.0% 


92.7% 


93.1% 


Width 


26.67 


26.25 


26.79 


27.77 


34.67 


Ordinary 
kriging 


PMSE 


144.58 


83.28 


57.16 


77.19 


218.03 


PMAE 


9.31 


6.89 


5.62 


6.51 


11.52 


Emp. cov. 


59.5% 


76.6% 


88.2% 


92.8% 


78.0% 


Width 


19.10 


19.73 


22.26 


30.26 


34.69 



Table 2 

Predictive Mean Square Error, Predictive Mean Absolute Error, empirical coverage and 
width of the 95% predictive interval for PM 2 . 5 for different values of the decay parameter 

under different models. 



Model 


Summary 
statistic 


Decay parameter 


0.125e-4 


0.125e-3 


0.125e-2 


0.125e-l 


0.125 


Independent 
downscaler 


PMSE 


13.46 


13.27 


12.95 


16.58 


36.43 


PMAE 


2.43 


2.42 


2.41 


2.65 


3.83 


Emp. cov. 


93.0% 


91.8% 


92.7% 


91.6% 


92.7% 


Width 


15.01 


14.29 


14.79 


15.05 


20.17 


Bivariate 
downscaler 


PMSE 


11.98 


11.78 


11.69 


14.43 


32.55 


PMAE 


2.33 


2.31 


2.29 


2.49 


3.55 


Emp. cov. 


92.7% 


93.6% 


94.0% 


93.2% 


93.5% 


Width 


14.02 


14.33 


14.53 


14.87 


18.63 


Kriging 


PMSE 


43.89 


21.99 


16.19 


25.65 


64.17 


PMAE 


4.37 


3.06 


2.54 


3.11 


5.37 


Emp. cov. 


52.2% 


72.5% 


89.7% 


94.8% 


89.9% 


Width 


7.32 


8.89 


14.06 


22.61 


25.73 
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Table 3 

Predictive Mean Square Error, Predictive Mean Absolute Error, empirical coverage, 
width of the 95% predictive interval, Continues Ranked Probability Score and Interval 
Score for both pollutants for the different models, kriging and cokriging at all the 65 

validation sites. 









Independent 


Bivariate 






Pollutant 




CMAQ 


downscaler 


downscaler 


Kriging 


Cokriging 




PMSE 


122.1 


55.9 


52.8 


57.2 


53.4 




PMAE 


8.5 


5.6 


5.4 


5.6 


5.4 


Ozone 


Emp. cov. 




92.9% 


93.0% 


88.2% 


94.3% 




Width 




27.3 


26.8 


22.3 


27.9 




CRPS 


8.5 


4.0 


3.9 


4.1 


5.2 




Interval 














score 




40.0 


38.9 


43.0 


37.1 




PMSE 


63.8 


13.0 


11.7 


16.2 


13.9 




PMAE 


5.1 


2.4 


2.3 


2.5 


2.4 


PM 2 .5 


Emp. cov. 




92.7% 


94.0% 


89.7% 


92.6% 




Width 




14.8 


14.6 


14.1 


12.8 




CRPS 


5.1 


1.8 


1.7 


1.9 


2.3 




Interval 














score 




19.2 


17.8 


21.1 


19.6 
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Table 4 

Predictive Mean Square Error, Predictive Mean Absolute Error, empirical coverage and 
width of the 95% predictive interval, Continuous Ranked Probability Score and Interval 

Score for both pollutants for the different models, kriging and cokriging at the 13 
validation sites that are at no more than 40 km of distance from the closest monitoring 

test site. 









Independent 


Bivariate 






Pollutant 




CMAQ 


downscaler 


downscaler 


Kriging 


Cokriging 




PMSE 


124.4 


57.4 


54.1 


57.9 


54.2 




PMAE 


8.5 


L 5.6 


5.4 


5.6 


5.5 


Ozone 


Emp. cov. 




92.7% 


92.7% 


88.4% 


94.3% 




Width 




27.3 


26.7 


22.2 


27.9 




CRPS 


8.5 


4.1 


4.0 


4.1 


5.2 




Interval 














score 




41.0 


39.9 


43.4 


37.7 




PMSE 


67.9 


13.3 


12.1 


16.9 


14.5 




PMAE 


5.2 


2.4 


2.3 


2.6 


2.4 


PM 2 . 5 


Emp. cov. 




92.8% 


94.4% 


89.9% 


92.9% 




Width 




13.8 


13.4 


12.1 


12.9 




CRPS 


5.2 


1.8 


1.7 


1.9 


2.3 




Interval 














score 




19.9 


18.3 


22.0 


18.8 
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Table 5 

Predictive Mean Square Error, Predictive Mean Absolute Error, empirical coverage and 
width of the 95% predictive interval, Continuous Ranked Probability Score and Interval 

Score for both pollutants for the different models, kriging and cokriging at the 52 
validation sites that are at more than 40 km of distance from the closest monitoring test 

site. 









Independent 


Bivariate 






Pollutant 




CMAQ 


downscaler 


downscaler 


Kriging 


Cokriging 




PMSE 


113.1 


50.1 


47.4 


54.4 


50.1 




PMAE 


8.2 


1 5 ' 3 


5.2 


5.6 


5.3 


Ozone 


Emp. cov. 




93.7% 


94.3% 


87.5% 


94.5% 




Width 




27.1 


26.5 


22.2 


27.8 




CRPS 


8.2 


3.9 


3.8 


4.1 


5.1 




Interval 














score 




36.1 


34.7 


41.4 


35.0 




PMSE 


48.5 


11.6 


10.0 


13.5 


11.4 




PMAE 


4.6 


2.4 


2.2 


2.4 


2.3 


PM 2 . 5 


Emp. cov. 




92.3% 


92.9% 


88.8% 


91.4% 




Width 




13.1 


12.8 


11.6 


12.3 




CRPS 


4.6 


1.7 


1.6 


1.8 


2.2 




Interval 














score 




16.6 


15.7 


17.6 


15.6 
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FIGURES 




Fig 1: Sites reporting concentration of ozone and PM2.5 used in our case 
study. Sites measuring only ozone are represented with dots, sites reporting 
only PM2.5 concentrations are represented with triangles, while sites mea- 
suring concentration for both pollutants are represented with squares. Black 
symbols are used to display sites used to fit the model, while grey symbols 
indicate validation sites. 
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Fig 2: Normal Q-Q plots of: (a) square root of observed ozone; (b) log of 
observed PM2.5. 
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FIGURES 




Fig 3: Time series of: (a) daily mean (open circles) and daily standard de- 
viation (black dots) of square root of observed ozone; (b) daily mean (open 
circles) and daily standard deviation (black dots) of log of observed PM2.5. 
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(a) Ozone - Intercept 



(b) PM2.S - Intercept 




(c) Ozone - Coefficient of CMAQ output for ozone 




Latitude 

30 32 34 36 38 40 




7 * .* • • :^ 


mm \ 




-90 -85 -80 -75 
Longitude 

(d) PM2.5 - Coefficient of CMAQ output for ozone 


l :: 




- 1.0 

- 0.5 oj 

a - 




- 0.0 


• 


- -0.5 





(e) Ozone - Coefficient of CMAQ output for PM2.5 



(f) PM2.5 - Coefficient of CMAQ output for PM2.5 





Fig 4: Spatial maps of the estimates of the coefficients of the linear regres- 
sions of the square root of observed concentration of ozone (panels (a), (c) 
and (e)) and of the logarithm of the observed concentration of PM2.5 (panels 
(b), (d) and (f)) on the square root of the CMAQ predicted concentration 
of ozone and on the log of the CMAQ predicted concentration of PM2.5: (a)- 
(b) intercept term; (c)-(d) coefficient of the CMAQ model output for ozone; 
(e)-(f) coefficient of the CMAQ model output for PM2.5. In all panels, the 
HMfiipffin n^bWri 5 caYA<ed Wem^^!fe 5] hc^filgd d re§p^e i ^d 2010 
the normalized covariates. 
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FIGURES 





Fig 5: (a) Observed ozone on June 25, 2002; (b) Predicted ozone as obtained 
from CMAQ on June 25, 2002; (c) Predicted ozone as obtained via kriging for 
June 25, 2002; (d) Predicted ozone as obtained from the bivariate downscaler 
model for June 25, 2002. 
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(a) 



(b) 





(c) 



(d) 



Fig 6: (a) Observed PM 2 . 5 on June 25, 2002; (b) Predicted PM 2 . 5 as ob- 
tained from CMAQ on June 25, 2002; (c) Predicted PM2.5 as obtained via 
kriging for June 25, 2002; (d) Predicted PM2.5 as obtained from the bivariate 
downscaler model for June 25, 2002. 
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FIGURES 
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Fig 8: Posterior predictive distribution for the difference in average ozone 
and PM2.5 concentration between Mid-Atlantic and Midwest, Mid-Atlantic 
and South and Midwest and South on June 25, 2002, as obtained using the 
bivariate downscaler. 



imsart-aoas ver. 2009/12/15 file: bidownscaler_paper_AOAS.tex date: April 8, 2010 



